
## Today's Packages and Data 🤗

```{r}
install.packages("GGally")
install.packages("devtools")
devtools::install_github("quinix45/FabioFun")
```

```{r}
# package 'devtools' must be installed first
library(GGally)
library(tidyverse)
library(FabioFun)
# I am setting a theme for all the plots
theme_set(theme_classic(base_size = 16, accent = "#7a0b80",
                        base_family = 'serif'))
```



```{r}
dat <- rio::import("https://fabio-setti.netlify.app/data/World_happiness_2024.csv")
```

## Predicting Country Happiness

```{r}
reg_free <- lm(Happiness_score ~ Freedom, dat)
summary(reg_free)
```


```{r}
reg_corr <- lm(Happiness_score ~ Corruption, dat)
summary(reg_corr)
```


## Well, That's just the correlation 🤷

```{r}
# here I am selecting the 3 columns that I wantr to plot
ggpairs(dat[,c(3, 7,9)])
```

## Multiple Independent Variables


```{r}
reg_multi <- lm(Happiness_score ~ Freedom + Corruption, dat)
summary(reg_multi)
```


## Individual regression plots


```{r}
ggplot(dat,  
 aes(x = Freedom, y = Happiness_score)) +
 geom_point() + 
    geom_smooth(method = "lm",  
                se = FALSE)
```

```{r}
ggplot(dat,  
 aes(x = Corruption, y = Happiness_score)) +
 geom_point() + 
    geom_smooth(method = "lm",  
                se = FALSE)
```


## Adding Multiple regression Lines

```{r}
ggplot(dat,  
 aes(x = Freedom, y = Happiness_score)) +
 geom_point() + 
    geom_smooth(method = "lm",  
                se = FALSE) +
    geom_abline(intercept = coef(reg_multi)[1], 
                slope = coef(reg_multi)[2])
```



```{r}
ggplot(dat,  
 aes(x = Corruption, y = Happiness_score)) +
 geom_point() + 
    geom_smooth(method = "lm",  
                se = FALSE)+
    geom_abline(intercept = coef(reg_multi)[1], 
                slope = coef(reg_multi)[3])
```

## A better Perspective


```{r}
library(FabioFun)

nice_3D_plot(y = dat$Happiness_score,
             x1 = dat$Freedom,
             x2 = dat$Corruption,
             dot_labels = dat$Country_name,
             axis_names = c("Happiness", 
                            "Freedom", 
                            "Corruption"))
```

## The Regression *Plane*


```{r}
nice_3D_plot(y = dat$Happiness_score,
             x1 = dat$Freedom,
             x2 = dat$Corruption,
             dot_labels = dat$Country_name,
             axis_names = c("Happiness", 
                            "Freedom", 
                            "Corruption"),
             reg_plane = TRUE)
```


## Added Variable Plots

```{r}
car::avPlots(reg_multi)
```


## Interpreting Multiple Regression

```{r}
summary(reg_multi)
```

## Checking Residuals


```{r}
car::residualPlots(reg_multi, fitted = FALSE, test = FALSE)
```



## The Scale of Our Variables and Interpretation


```{r}
coef(reg_multi)
```


```{r}
range(dat$Happiness_score)
range(dat$Freedom)
range(dat$Corruption)
```


```{r}
dat$Happiness_std <- scale(dat$Happiness_score)[,1]
dat$Freedom_std <- scale(dat$Freedom)[,1]
dat$Corruption_std <- scale(dat$Corruption)[,1]
```

## Standardized Multiple Regression

```{r}
reg_std <- lm(Happiness_std ~ Freedom_std + Corruption_std, dat)
summary(reg_std)
```

## APA Style Table for Regression Results


```{r}
library(FabioFun)
reg_tab <- regression_flextable(reg_std, 
                                var_names = c("Freedom", "Corruption"), 
                                intercept = FALSE)
reg_tab

# save to Word document (will be saved to current working directory)

# uncomment to save 
# flextable::save_as_docx(reg_tab, path = "table.docx")
```



## Categorical Variables 

```{r}
table(dat$Region)
```


## Dummy Coding With Factor Variables


```{r}
dat_filt <- dat %>% 
   filter(Region %in% c("Africa", "Asia", "Western Europe", "Eastern Europe"))
```


```{r}
dat_filt$Region <- factor(dat_filt$Region)
```


```{r}
contrasts(dat_filt$Region)
```

## Regression With a Categorical Variable

```{r}
reg_cat <- lm(Happiness_score ~ Region, dat_filt)
summary(reg_cat)
```


```{r}
dat_filt %>% 
   group_by(Region) %>% 
    summarise(Happ_mean = mean(Happiness_score))
```

## The Graphical intuition


```{r}
dat_coef <- data.frame(c("Intercept", "Asia", "Eastern Europe", "Western Europe"), 
                        as.numeric(round(coef(reg_cat), 2)))
colnames(dat_coef) <- c("Coefficient", "Value")
dat_coef

```


```{r}

happ_means <- dat_filt %>% 
   group_by(Region) %>% 
    summarise(Happ_mean = mean(Happiness_score))

# reference mean 
ref <- as.numeric(happ_means[1,2])

ggplot(happ_means, aes(x = Region, y = Happ_mean)) +
  geom_point(size = 4, color = "#7a0b80") +
  geom_text(aes(label = round(Happ_mean, 2)), vjust = -1) +
  geom_hline(yintercept = ref, lty = 3) +
  geom_segment(x = "Asia", 
               xend = "Asia", 
               y = ref,
               yend = as.numeric(happ_means[2,2]), lty = 3) +
 geom_segment(x = "Eastern Europe", 
               xend = "Eastern Europe", 
               y = ref,
               yend = as.numeric(happ_means[3,2]), lty = 3) +
geom_segment(x = "Western Europe", 
               xend = "Western Europe", 
               y = ref,
               yend = as.numeric(happ_means[4,2]), lty = 3) +
annotate("text", x = 1,  y = 4.2, label = "Reference Point") +
annotate("text", x = 2.2,  y = 5.3, label = paste("+", round(coef(reg_cat)[2], 2))) +
annotate("text", x = 3.2,  y = 5.3, label = paste("+", round(coef(reg_cat)[3], 2))) +
annotate("text", x = 4.2,  y = 5.3, label = paste("+", round(coef(reg_cat)[4], 2))) +
ylim(c(4, 6.9)) +
ylab("Mean Region Happiness")
```

## Why Does this work?

```{r}
intercept <- coef(reg_cat)[1]
slopes <- coef(reg_cat)[-1]

# the intercept + all slopes multiplied by 0. Will give the mean happiness for Africa
intercept + sum(slopes*c(0,0,0))
```


```{r}
# mean happiness of western Europe (equivalent to doing 4.40 + 2.35)
intercept + sum(slopes*c(0,0,1))
```

## The Moral of the Story


```{r}
contrasts(dat_filt$Region)[contrasts(dat_filt$Region) == 1] <- 3
contrasts(dat_filt$Region)
```


```{r}
reg_cat <- lm(Happiness_score ~ Region, dat_filt)
summary(reg_cat)
```